occurrence <- function (component) { nrow (component) }

correlation <- function (c1, c2) { nrow (c1) / nrow (c2) }

distribution <- function (component)
{
    class_counts <- component %>% add_count (class) %>% extract2 ("n")
    coords <- st_coordinates (component)
    idists <- coords %>% dist %>% as.matrix %>% divide_by (1)
    diag (idists) <- 0
    m_i <- Moran.I (class_counts, idists) %>% extract2 ("observed")

    return (m_i)
}

proximity <- function (c1, c2)
{
    if (nrow (c2) == 0)
    {
        dists <- NA
    } else
    {
        c1_pnt <- st_centroid (c1)
        dists <- st_distance (c1_pnt, c2) %>% as.numeric %>% mean
    }

    return (dists)
}

mean_dist_mat <- function (c1)
{
    d <- NA
    if (nrow (c1) > 1)
        d <- st_distance (c1) %>% as.numeric %>% mean
    d
}
tc <- function (func)
{
    res <- try (func, silent = TRUE)
    if (inherits (res, "try-error"))
        res <- NA

    return (res)
}

cnames_res <- c ("occurrence_shops_500m", "occurrence_amenities_500m",
                  "occurrence_shops", "occurrence_amenities",
                  "occurrence_junctions_5000m", "occurrence_highways_5000m",
                  "distribution_shop", "distribution_amenities",
                  "distribution_shop_500ms", "distribution_amenities_500m",
                  "proximity_shops_500m", "proximity_amenities_500m",
                  "proximity_shops", "proximity_amenities",
                  "proximity_transport_500m", "proximity_junction_5km",
                  "proximity_highway_5km", "correlation_shops_amenities")
res <- matrix (NA, nrow = nrow (malls), ncol = length (cnames_res)) %>%
    as.data.frame
names (res) <- cnames_res

if (!results_gpkg %in% dir (dir_out))
{
    for (i in seq_len (nrow (malls)))
    {
        print (i)
        mall <- malls %>% extract (i, ) %>% st_buffer (5)
        mall_500m <- malls_500m %>% extract (i, ) %>% extract ("geom")
        mall_5000m <- malls_5000m %>% extract (i, ) %>% extract ("geom")

        amenities_500m_i <- clip_pts_by_poly (amenities, mall_500m)
        amenities_i <- clip_pts_by_poly (amenities, mall)
        highways_i <- clip_pts_by_poly (highways, mall_5000m)
        junctions_i <- clip_pts_by_poly (junctions, mall_5000m)
        public_transport_i <- clip_pts_by_poly (public_transport, mall_500m)
        shops_500m_i <- clip_pts_by_poly (shops, mall_500m)
        shops_i <- clip_pts_by_poly (shops, mall)

        res [i, 1] <- tc (occurrence (shops_500m_i))
        res [i, 2] <- tc (occurrence (amenities_500m_i))
        res [i, 3] <- tc (occurrence (shops_i))
        res [i, 4] <- tc (occurrence (amenities_i))
        res [i, 5] <- tc (occurrence (junctions_i))
        res [i, 6] <- tc (occurrence (highways_i))
        #res [i, 7] <- tc (distribution (shops_i))
        #res [i, 8] <- tc (distribution (amenities_i))
        #res [i, 9] <- tc (distribution (shops_500m_i))
        #res [i, 10] <- tc (distribution (amenities_500m_i))
        res [i, 11] <- tc (proximity (mall, shops_500m_i))
        res [i, 12] <- tc (proximity (mall, amenities_500m_i))
        res [i, 13] <- tc (proximity (mall, shops_i))
        res [i, 14] <- tc (proximity (mall, amenities_i))
        res [i, 15] <- tc (proximity (mall, public_transport_i))
        res [i, 16] <- tc (proximity (mall, junctions_i))
        res [i, 17] <- tc (proximity (mall, highways_i))
        res [i, 18] <- tc (correlation (shops_i, amenities_i))
    }
    res$correlation_shops_amenities [(res [, 18] == Inf)] <- NA
    malls_result <- bind_cols (malls, res)
    res_min <- sapply (res, min, na.rm = TRUE) %>% t %>% unname
    res_mean <- sapply (res, mean, na.rm = TRUE) %>% t %>% unname
    res_median <- sapply (res, median, na.rm = TRUE) %>% t %>% unname
    res_max <- sapply (res, max, na.rm = TRUE) %>% t %>% unname
    res_sd <- sapply (res, sd, na.rm = TRUE) %>% t %>% unname
    res_cv <- mod ((100 * res_sd / res_mean), 100)

    

    res_stats <- rbind (res_min, res_mean, res_median, res_max, res_sd, res_cv)
    colnames (res_stats) <- cnames_res
    rownames (res_stats) <- c ("min", "mean", "median", "max", "sd", "cv")

    write.table (res, paste0 (dir_out_csv, "results.csv"), row.names = FALSE,
                 sep = ",", quote = FALSE)
    write.csv (res_stats, paste0 (dir_out_csv, "results_stats.csv"),
               quote = FALSE)

    st_write (malls_result, paste0 (dir_out, results_gpkg), quiet = TRUE)
}
malls_result <- st_read (paste0 (dir_out, results_gpkg), quiet = TRUE)
osm_l <- osmdata::getbb ("London", format_out = "polygon") %>% extract2 (1) %>%
    extract2 (1)
wkt <- paste0 ("SRID=4326;POLYGON((",
                              paste (paste (osm_l [, 1], osm_l [, 2]), collapse
                                     = ","),
                              "))")
aoi <- st_as_sfc (wkt) %>% st_set_crs (4326) %>% st_transform (crs_utm_london)
aoi_sp <- as (aoi, "Spatial") %>% geometry
pts <- spsample (aoi_sp, type = "hexagonal", cellsize = 500, offset = c (0, 0))
hexgrid <- HexPoints2SpatialPolygons (pts) %>% as ("sf")

cnames_hex <- c ("occurrence_shops", "occurrence_amenities",
                 "occurence_bus_stops", "dmat_shops",
                 "dmat_amenities", "dmat_shops_amenities",
                 "correlation_shops_amenities")

res_hex <- matrix (NA, nrow = nrow (hexgrid), ncol = length (cnames_hex)) %>%
    as.data.frame
names (res_hex) <- cnames_hex
shops %<>% clip_pts_by_poly (aoi)
amenities %<>% clip_pts_by_poly (aoi)

int_hex_public_t <- st_intersects (hexgrid, public_transport)
int_hex_shop <- st_intersects (hexgrid, shops)
int_hex_am <- st_intersects (hexgrid, amenities)

if (!results_hex_gpkg %in% dir (dir_out))
{
    for (i in seq_len (nrow (hexgrid)))
    {
        print (paste0 (i, "/", nrow(hexgrid)))

        shop_ids <- int_hex_shop %>% extract2 (i)
        amenity_ids <- int_hex_am %>% extract2 (i)
        pt_ids <- int_hex_public_t %>% extract2 (i)

        shops_i <- shops %>% extract (shop_ids, )
        public_transport_i <- public_transport %>% extract (pt_ids, )
        amenities_i <- amenities %>% extract (amenity_ids, )

        occ_sh <- tc (occurrence (shops_i))
        occ_am <- tc (occurrence (amenities_i))
        res_hex [i, 1] <- occ_sh
        res_hex [i, 2] <- occ_am
        res_hex [i, 3] <- tc (occurrence (public_transport_i))
        res_hex [i, 4] <- tc (mean_dist_mat (shops_i))
        res_hex [i, 5] <- tc (mean_dist_mat (amenities_i))
        res_hex [i, 6] <- tc (mean_dist_mat (rbind (shops_i, amenities_i)))
        res_hex [i, 7] <- occ_sh / occ_am
    }
    res_hex[, 7] [is.infinite (res_hex [, 7])] <- NA

    hexgrid <- bind_cols (hexgrid, res_hex)
    res_hex_min <- sapply (res_hex, min, na.rm = TRUE) %>% t %>% unname
    res_hex_mean <- sapply (res_hex, mean, na.rm = TRUE) %>% t %>% unname
    res_hex_median <- sapply (res_hex, median, na.rm = TRUE) %>% t %>% unname
    res_hex_max <- sapply (res_hex, max, na.rm = TRUE) %>% t %>% unname
    res_hex_sd <- sapply (res_hex, sd, na.rm = TRUE) %>% t %>% unname
    res_hex_cv <- mod ((100 * res_hex_sd / res_hex_mean), 100)

    res_hex_stats <- rbind (res_hex_min, res_hex_mean, res_hex_median,
                            res_hex_max, res_hex_sd, res_hex_cv)
    colnames (res_hex_stats) <- cnames_hex
    rownames (res_hex_stats) <- c ("min", "mean", "median", "max", "sd", "cv")

    write.table (res_hex, paste0 (dir_out_csv, "results_hex.csv"),
                 row.names = FALSE, sep = ",", quote = FALSE)
    write.csv (res_hex_stats, paste0 (dir_out_csv, "results_hex_stats.csv"),
               quote = FALSE)

    st_write (hexgrid, paste0 (dir_out, results_hex_gpkg), quiet = TRUE)
}
hexgrid <- st_read (paste0 (dir_out, results_hex_gpkg), quiet = TRUE)
hexgrid$score <- 0
hexgrid$mall_contained <- 0
fs <- hexgrid$occurrence_shops >= 4

p_am <- fs & (hexgrid$occurrence_amenities > 0)
p_d_shop <- fs & (hexgrid$dmat_shops >= 5.8 )
p_d_am <- fs & (hexgrid$dmat_amenities >= 6.5)
p_corr <- fs & (hexgrid$correlation_shops_amenities > 0.6)
p_am [is.na (p_am)] <- FALSE
p_d_shop [is.na (p_d_shop)] <- FALSE
p_d_am [is.na (p_d_am)] <- FALSE
p_corr [is.na (p_corr)] <- FALSE

hexgrid$score <- colSums (rbind (fs, p_am, p_d_shop, p_d_am, p_corr))
malls$score <- 0

for (i in seq_len (nrow (malls)))
{
    mall <- malls [i, ]
    hex_idx <- st_intersects (mall, hexgrid) %>% unlist
    hexgrid$mall_contained [hex_idx] <- hex_idx
    hex_tile <- hexgrid [hex_idx, ]
    malls$score [i] <- hex_tile %>% extract2 ("score") %>% mean
}
mean (malls$score, na.rm = TRUE)

[1] 4.488782

mean (hexgrid$score, na.rm = TRUE)

[1] 0.8563491

for (cname in cnames_hex)
    plot (hexgrid[, cname])

## Warning in classInt::classIntervals(na.omit(values), min(nbreaks, n.unq), :
## var has infinite values, omitted in finding classes

plot (hexgrid[, "score"])